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ABSTRACT 

The  computational  complexity  of  track  initiation,  also  known  as  initial  orbit  determination  or  IOD,  using  only  angle 
measurements  is  polynomial  in  the  number  of  observations.  However,  the  polynomial  degree  can  be  high,  always  at 
least  cubic  and  commonly  quartic  or  higher.  Therefore,  practical  implementations  require  attention  to  the  scalability 
of  the  algorithms,  when  one  is  dealing  with  the  very  large  number  of  observations  from  large  surveillance 
telescopes.  We  address  two  broad  categories  of  algorithms.  The  first  category  includes  and  extends  the  classical 
methods  of  Laplace  and  Gauss,  as  well  as  the  more  modern  method  of  Gooding,  in  which  one  solves  explicitly  for 
the  apparent  range  to  the  target  in  terms  of  the  given  data.  We  find  that  the  orbit  solutions  (data  association 
hypotheses)  can  be  ranked  by  means  of  a  concept  we  call  persistence ,  in  which  a  simple  statistical  measure  of 
likelihood  is  based  on  the  frequency  of  occurrence  of  combinations  of  observations  in  consistent  orbit  solutions. 
However,  range-solution  methods  can  be  expected  to  perform  poorly  if  the  initial  orbit  solutions  of  most  interest  are 
not  well  conditioned.  The  second  category  of  algorithms  addresses  this  difficulty.  Instead  of  solving  for  range,  these 
methods  attach  a  set  of  range  hypotheses  to  each  measured  line  of  sight.  Then  all  pair-wise  combinations  of 
observations  are  considered  and  the  family  of  Lambert  problems  is  solved  for  each  pair.  These  algorithms  also  have 
polynomial  complexity,  though  now  the  complexity  is  quadratic  in  the  number  of  observations  and  also  quadratic  in 
the  number  of  range  hypotheses.  We  offer  a  novel  type  of  admissible-region  analysis,  constructing  partitions  of  the 
orbital  element  space  and  deriving  rigorous  upper  and  lower  bounds  on  the  possible  values  of  the  range  for  each 
partition.  This  analysis  allows  us  to  parallelize  with  respect  to  the  element  partitions  and  to  reduce  the  number  of 
range  hypotheses  that  have  to  be  considered  in  each  processor  simply  by  making  the  partitions  smaller.  We  present 
numerical  results  based  on  simulated  data  sets. 


1.  MOTIVATION 

The  advent  of  high-sensitivity,  high-capacity  optical  sensors  for  space  surveillance  presents  us  with  interesting  and 
challenging  tracking  problems.  Accounting  for  the  origin  of  every  detection  made  by  such  systems  is  generally 
agreed  to  belong  to  the  “most  difficult”  category  of  tracking  problems.  Especially  in  the  early  phases  of  the  tracking 
scenario,  when  a  catalog  of  targets  is  being  compiled,  or  when  many  new  objects  appear  in  space  because  of  on-orbit 
explosion  or  collision,  one  faces  a  combinatorially  large  number  of  orbit  (data  association)  hypotheses  to  evaluate. 
The  number  of  hypotheses  is  reduced  to  a  more  feasible  number  if  observations  close  together  in  time  can,  with  high 
confidence,  be  associated  by  the  sensor  into  extended  tracks  on  single  objects.  Most  current  space  surveillance 
techniques  are  predicated  on  the  sensor  systems’  ability  to  form  such  tracks  reliably.  However,  the  required 
operational  tempo  of  space  surveillance,  the  very  large  number  of  objects  in  Earth  orbit  and  the  difficulties  of 
detecting  dim,  fast-moving  targets  at  long  ranges  means  that  individual  sensor  track  reports  are  often  inadequate  for 
computing  initial  orbit  hypotheses.  In  fact,  this  situation  can  occur  with  optical  sensors  even  when  the  probability  of 
detection  is  high.  For  example,  the  arc  of  orbit  that  has  been  observed  may  be  too  short  or  may  have  been  sampled 
too  sparsely  to  allow  well-conditioned,  usable  orbit  estimates  from  single  tracks.  In  that  case,  one  has  no  choice  but 
to  solve  a  data  association  problem  involving  an  unknown  number  of  targets  and  many  widely  spaced  observations 
of  uncertain  origin.  In  the  present  paper,  we  are  motivated  by  this  more  difficult  aspect  of  the  satellite  cataloging 
problem.  However,  the  results  of  this  analysis  may  find  use  in  a  variety  of  less  stressing  tracking  applications. 

One  method  for  track  initiation  from  the  asteroid  community  is  the  use  of  admissible  regions  [1]  [2].  This  work  has 
been  adapted  for  use  in  tracking  Earth  orbiting  objects  by  Dan  Scheeres  and  his  group  at  University  Colorado  - 
Boulder  [3].  However,  this  methodology  was  developed  by  astronomers  who  had  access  to  many  nights  of  data  for 
relatively  slow  moving  objects  (from  an  Earth  point  of  view),  and  it  requires  simultaneous  angle  and  angle  rate  data 
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or  the  ability  to  synthesize  those  data  from  observed  quantities.  For  GEO  objects,  we  typically  have  only  angles 
data.  Our  ability  to  estimate  angle  rate  is  limited.  Hence,  the  usefulness  of  this  method  for  our  purpose  is  equally 
limited  but  nevertheless  has  shown  some  promise  in  analyses  by  Scheeres’  group.  Without  having  angle  rate 
available,  we  still  have  the  option  to  use  the  angles  data  in  conjunction  with  guesses  for  range.  The  range  guesses 
could  be  given  in  a  predetermined  sequence  over  some  allowable  set  of  values  or  could  be  sampled  randomly.  In 
either  case,  the  allowable  set  of  range  values  would  be  bounded  by  constraints  derived  from  the  allowable  set  of 
orbital  elements  used  for  the  track  initiation.  Then  Lambert’s  Theorem  would  be  used  to  generate  a  set  of  initial 
orbit  hypotheses  for  each  pair  of  measured  lines  of  sight.  Lambert’s  Theorem  states  that  the  time  of  flight  between 
two  positions  in  space  depends  only  on  the  semi-major  axis  of  the  orbit  and  the  position  vectors  themselves.  Hence, 
two  positions  and  the  time  of  flight  between  them  completely  specify  the  orbit.  Alternatively,  we  can  seek  to  solve 
for  the  range  at  each  observation  time  by  imposing  dynamical  constraints.  This  approach  has  higher  computational 
complexity  but  yields  higher  confidence  earlier  in  the  orbit  determination  process.  To  speed  up  this  process,  we  seek 
to  reduce  the  number  of  hypotheses  a  priori.  First,  we  apply  geometrical  boundaries  on  range  hypotheses.  These 
boundaries  may  be  somewhat  overbroad  to  ensure  that  the  actual  orbit  lies  within  the  search  space.  Second,  we  solve 
iteratively  for  ranges  and  use  the  geometrical  boundaries  to  screen  out  bad  solutions. 

How  should  we  limit  the  number  of  range  hypotheses  to  make  the  total  number  of  candidate  orbits  manageable 
while  also  generating  candidates  that  are  likely  to  correspond  to  real  orbits  of  interest?  Let  us  seek  to  generate 
hypotheses  for  orbits  that  lie  only  in  a  bounded  region  of  semimajor  axis  a,  eccentricity  e,  inclination  /  and  right 
ascension  of  the  ascending  node  H,  namely,  within  a  partition  specified  by  the  intervals  [aMIN,aMAX], 
[^min  t  ^maxL  [?min  ^max]  and  [HMin  >^maxL  F°r  the  purposes  of  this  discussion,  we  leave  the  other  orbital 
elements  unconstrained.  It  will  turn  out  that  these  elements  constrain  the  possible  values  of  range  in  simple  ways 
without  our  having  any  recourse  to  angle  rate  information.  Then,  to  the  extent  that  we  can  restrict  the  generation  of 
hypothetical  orbits  to  a  specified  partition  of  the  space  of  orbital  elements,  we  have  parallelized  the  task  of  building 
a  catalog  of  objects  detected  within  that  partition.  The  reason  is  that  any  partition  of  the  space  of  orbital  elements, 
including  the  whole  space  itself,  can  be  sub-divided  into  smaller  partitions,  and  each  sub-partition  can  be  handled 
independently.  In  the  approach  outlined  here,  all  the  observations  would  have  to  be  considered  for  each  sub-partition 
of  the  space  of  orbit  elements.  However,  by  constructing  upper  and  lower  bounds  on  range  for  each  sub-partition  of 
the  element  space,  we  limit  the  number  of  range  hypotheses  that  have  to  be  considered  for  each  sub-partition.  This 
approach  allows  us  to  consider  a  manageable  number  of  range  hypotheses  for  each  sub-partition  before  we  generate 
candidate  orbits,  simply  by  making  the  sub-partitions  small  enough,  so  that  the  overall  computation  is  feasible. 

2.  RANGE  BOUNDARY  FORMULATION 

Assume  that  we  have  a  pair  of  line-of-sight  unit 
vectors  ux  and  u2  ,  measured  at  time  tx  at  station 
position  Rx  and  time  t2  at  station  position  R2, 
respectively.  Assume,  without  loss  of  generality, 
that  t2  >  G.  We  want  to  form  the  hypothesis  that 
these  two  observations  are  associated  with  the 
same  space  object.  To  this  end,  we  attach  a  set  of 
hypothetical  range  values,  =  1,2, ...)  and 

{P2,n  >n  =  1,2,...}  respectively,  to  each  of  these 
measured  unit  vectors  and  then  generate  candidate 
orbits  by  solving  Lambert’s  problem  for  each  of 
the  pair-wise  combinations  of  hypothetical  orbital 
position  vectors  rlm  and  r2  n. 

We  are  seeking  explicit  bounds  on  range  that  can 
be  applied  for  each  individual  angle-based 
observation,  or  at  most  to  pairs  of  angle-based 
observations.  Even  with  the  further  restriction  that 
hypothetical  orbits  be  elliptical  and  Keplerian  (which  we  accept),  it  may  not  be  obvious  that  efficient  bounds  having 
these  properties  can  be  obtained.  Exact  bounds  would  have  to  be  based  on  some  admissible-region  analysis  of  the 


Figure  1:  Vector  Triangles 


type  developed  by  Milani,  Tommei,  Scheeres,  Maruskin,  Fujimoto  and  others.  For  example,  denoting  the 
gravitational  parameter  by  fi  ,  we  write  the  first  integrals  of  Keplerian  motion  as 

energy:  E  =  (r  ■  r)/2  -  ju/||r||  =  -  n/2a  (i) 

angular  momentum:  h  =  r  x  r  (ii) 

Laplace  vector:  pe  =  r  x  (r  x  r)  —  jur/||r||  (iii) 

These  can  be  evaluated  with  the  vector  triangle  relation  r  =  R  +  pu  and  its  time  derivative  r  =  R  +  pu  +  pu  for 
each  observation.  Then,  for  each  observation,  we  can  define  admissible  regions  in  the  (p,  p)  plane  for  each  partition 
in  the  space  of  elements  by  means  of  inequalities  such  as 

—  p/(2aMIN)  <E<  —  p/(2aMAx)  (iy) 

cos /MAX  <  (h/||h||)  ■  k  <  cos /MIN  (v) 

eMlN  ^  llell  ^  eMAX  (v0 

Here  k  is  the  north  polar  unit  vector  in  the  Earth-centered  inertial  frame.  For  each  observation,  the  values  of  range 
and  range  rate  that  satisfy  these  inequalities  will  result  in  orbits  that  lie  only  within  the  given  partition  of  the  space  of 
elements.  DeMars  and  Jah  [6]  have  shown  what  the  admissible  regions  look  like  for  partitions  of  semimajor  axis  and 
eccentricity  by  a  numerical  treatment  of  the  above  inequalities.  Maruskin,  et  al .,  [7]  have  shown  how  the  admissible 
regions  evolve  in  time  and  how  the  overlap  of  the  admissible  regions  for  different  observations  can  help  solve  the 
data  association  problem.  However,  even  though  expressions  (i)  through  (vi)  can  be  reduced  to  polynomial  forms  in 
range  and  range  rate,  each  relation  is  coupled  in  both  variables  and  the  polynomial  degree  is  high,  preventing  us 
from  obtaining  explicit  expressions  for  range  and  range  rate  in  terms  of  the  given  data.  Moreover,  the  usual 
admissible-region  analysis  leads  nowhere  if  angle  rates  are  not  available.  For  example,  the  track-initiation  method  of 
DeMars,  et  al .,  [8],  involving  multiple  hypotheses  on  range  and  range  rate,  requires  both  angle  and  angle  rate  values. 

In  the  present  analysis,  we  take  a  geometric  and  kinematic  approach  that  does  lead  to  explicit  upper  and  lower 
bounds  on  the  possible  values  of  range  for  each  observation  or  pair  of  observations,  given  only  angle  data  at  discrete 
times.  In  fact,  we  find  several  inequalities  that  must  be  satisfied  simultaneously,  and  we  can  take  the  most  restrictive 
superposition  of  the  different  bounds  as  our  working  result.  In  case  angle  rates  are  available,  we  can  obtain  explicit 
upper  and  lower  bounds  on  range  rate,  as  well  as  additional  bounds  on  range.  It  may  happen  that,  for  a  given 
observation,  there  are  no  values  of  the  range  or  range  rate  that  lead  to  orbits  within  the  given  element-space 
partition,  so  that  the  observation  can  be  eliminated  from  further  consideration.  We  obtain  explicit  conditions  for  the 
existence  of  possible  values  of  range  and  range  rate,  in  terms  of  the  observation  itself. 

The  price  for  obtaining  explicit  bounds  on  range  and  range  rate  is  that  the  bounds  are  not  exact  but  somewhat 
conservative.  Although  every  orbit  within  the  element-space  partition  corresponds  to  values  of  range  and  range  rate 
that  lie  within  the  bounds  given  here,  some  values  of  range  or  range  rate  that  satisfy  the  bounds  may  lead  to  orbits 
that  lie  outside  the  given  partition.  This  situation  represents  inefficiency  in  the  parallelization  of  building  the  catalog: 
nearly  the  same  candidate  orbits  near  the  boundaries  of  the  element-space  partitions  may  be  generated  in  both  of  the 
adjacent  partitions,  if  the  range  or  range  rate  hypotheses  are  planted  densely  enough.  On  the  other  hand,  no 
candidate  orbits  within  the  given  element-space  partitions  will  be  missed  because  of  the  bounds  given  here.  The 
extent  and  cost  of  the  inefficient  duplication  of  candidate  orbits  will  depend  on  the  particular  datasets  and  element 
partitions  of  interest,  and  may  require  further  study.  In  practice,  of  course,  within  any  element  partition,  any  of  these 
extra  orbit  hypotheses  can  be  either  kept  or  discarded.  If  they  are  kept,  one  would  have,  at  most,  a  bookkeeping 
problem  of  transferring  the  extra  orbits  to  the  correct  element  partition.  The  trade-off  in  this  case  is  that  merely 
moving  data  between  processors  always  takes  time. 

First,  we  present  bounds  on  range  that  must  hold  for  each  observed  line  of  sight.  Assuming  that  all  orbits  of  interest 
are  elliptical,  require  that  the  orbital  radii  lie  between  the  maximum  specified  apogee  and  the  minimum  specified 
perigee: 

aMIN  (1  —  eMAx)  ^  llrll  ^  aMAx(l  +  eMAx)  (1) 

The  values  of  range  that  correspond  to  these  limits  on  orbital  radius  can  be  found  explicitly  using  the  vector  triangle 
relationship  r  =  R  +  pu  .  Squaring  terms  to  remove  the  radical,  we  have 


“minC1  -  eMAx)2  <  R  ■  R  +  2(R  ■  u)p  +  p2  <  a£,AX(l  +  eMAX)2 


(2) 


Consider  the  perigee  and  apogee  cases  separately.  For  the  perigee  case,  we  require  the  orbital  radius  to  be  no  smaller 
than  the  smallest  allowable  perigee  radius: 

P2  +  2(R  ■  u)p  -  K,n  (1  -  eMAX)2  -  R  ■  R]  >  0  (3) 

We  will  have  real  roots  if  and  only  if  the  argument  of  the  square  root  of  the  solution  is  non-negative: 

aMiisf(l  —  ^max)2  >  R  ■  [R  —  (R  ■  u)u]  (4) 

If  no  real  roots  of  the  quadratic  expression  (3)  exist,  then  we  can  immediately  discard  the  current  observation  and 
form  no  hypotheses  with  it.  The  reason  is  that  no  value  of  the  range  will  be  found  for  this  observation,  which  is 
consistent  with  the  specified  intervals  of  the  orbital  elements. 

Descartes’  rule  of  signs  tells  us  the  number  of  positive  real  roots.  If  the  third  coefficient  in  the  quadratic  form  (3)  is 
negative,  that  is,  if  IN(1  —  eMAX)2  >  R  ■  R  ,  then,  regardless  of  the  sign  of  the  second  coefficient  2(R  ■  u),  we 
will  have  one  positive  real  root  and  necessarily  also  one  negative  root.  Because  the  quadratic  is  concave-up,  the 
inequality  is  satisfied  to  the  left  of  the  negative  root  and  to  the  right  of  the  positive  root.  We  can  ignore  the  negative 
root  and  all  values  to  the  left  of  it,  because  we  require  a  priori  that  range  values  to  be  non-negative.  What  remains  is 
a  positive  lower  limit  on  the  possible  values  of  range: 

p  >  -(R  ■  u)  + 

It  is  worth  noting  that,  for  Earth-bound  stations,  the  third  coefficient  of  (3)  will  essentially  always  be  negative 
because  the  inequality  IN(1  —  eMAX)2  >  R  ■  R  is  approximately  the  condition  that  the  minimum  allowable 
perigee  radius  be  larger  than  the  Earth  radius.  Moreover,  the  second  coefficient  2(R  ■  u)  will  essentially  always  be 
positive  because  observations  have  to  be  taken  above  the  local  horizontal  plane  at  some  positive  local  elevation 
angle. 

For  space-based  observing  stations,  it  is  possible  that  neither  of  these  circumstances  would  be  true:  the  station’s 
orbital  position  may  be  higher  than  the  minimum  specified  perigee  radius,  or  observations  may  be  taken  at  negative 
local  elevation  angles,  or  both.  If  the  third  coefficient  in  (3)  is  positive,  that  is,  if  IN(1  —  eMAX)2  <  R  ■  R  ,  then 
the  quadratic  will  have  either  no  positive  real  roots  or  two  positive  real  roots,  depending  on  the  sign  of  the  second 
coefficient.  If,  furthermore,  the  second  coefficient  in  (3)  is  positive,  that  is,  if  (R  ■  u)  >  0  ,  then  we  have  no  positive 
real  roots,  but  only  a  pair  of  negative  roots.  Because  the  quadratic  is  concave-up,  the  inequality  (3)  is  satisfied  to  the 
left  of  the  more  negative  root  and  to  the  right  of  the  less  negative  root.  However,  since  we  require  a  priori  that  range 
values  be  non-negative,  we  are  left  merely  with  the  condition  that  p  >  0  .  If  the  second  coefficient  is  negative,  that 
is,  (Ru)  <  0  ,  meaning  that  the  observation  is  taken  at  negative  local  elevation  angle,  then  the  quadratic  will  have 
two  positive  real  roots.  Because  the  quadratic  is  concave -up,  the  inequality  (3)  will  be  satisfied  to  the  left  of  the 
smaller  root,  that  is,  between  p  =  0  and  the  smaller  root,  and  also  to  the  right  of  the  larger  root.  In  this  case,  we  have 
two  disjoint  intervals  of  range,  one  finite  and  one  semi-infinite,  over  which  range  hypotheses  will  satisfy  the  perigee 
constraint: 

0  <  p  <  -(R  ■  u)  - 


p  >  —  (R  ■  u)  + 

Now  we  consider  the  apogee  case  and  seek  to  derive  results  that  are  analogous  to  those  above.  The  apogee  case  will 
provide  us  with  conditions  on  values  of  the  range  that  are  complementary  to  those  of  the  perigee  case.  Since  both 
sets  of  conditions  must  be  satisfied  simultaneously,  we  can  take  the  most  restrictive  superposition  of  all  conditions 
on  range  to  define  the  set  of  values  over  which  we  must  form  range  hypotheses. 


For  the  apogee  case,  we  have  from  the  inequality  (2)  that  the  orbital  radius  must  be  no  larger  than  the  maximum 
allowable  apogee  radius: 


P2  +  2(R  ■  u)p  —  [ay[ ax(1  +  ^max)2  —  R  '  R]  <  0  (8) 

We  will  have  real  roots  if  and  only  if  the  argument  of  the  square  root  in  the  solution  is  non-negative: 

aMAx(l  +  £Max)2  >  R  ■  [R  —  (R  ■  u)u]  (9) 

If  no  real  roots  exist,  then  we  can  immediately  discard  the  observation  and  form  no  hypotheses  with  it.  The  reason  is 
that  no  value  of  the  range  will  be  found  for  this  observation,  which  is  also  consistent  with  the  specified  intervals  of 
the  orbital  elements.  Assuming  that  we  have  real  roots  in  equation  (8),  we  use  Descartes’  rule  of  signs  to  determine 
the  number  of  positive  real  roots.  If  the  third  coefficient  in  the  quadratic  form  (8)  is  negative,  that  is,  if  a^Ax(l  + 
^max)2  >  R  '  R,  then,  regardless  of  the  sign  of  the  second  coefficient  2(R  ■  u)  ,  we  will  have  one  positive  real  root 
and  necessarily  also  one  negative  root.  Because  the  quadratic  is  concave -up,  the  inequality  (8)  is  satisfied  between 
the  roots.  Moreover,  we  require  a  priori  that  range  values  be  non-negative,  so  we  can  say  without  loss  of  generality 
that  the  inequality  will  be  satisfied  between  p  =  0  and  the  positive  real  root.  The  result  is  that  we  have  an  upper 
bound  on  the  possible  values  of  range: 


0  <  p  <  -(R  ■  u)  + 


(R  ■  u)2  +  [a 


2 

MAX 


(1  +  eMAx)2  —  R  '  R] 


(10) 


As  before,  for  Earth-bound  stations,  the  third  coefficient  will  essentially  always  be  negative  because  the  inequality 
aMAx(l  +  ^max)2  >  R  '  R  is  approximately  the  condition  that  the  maximum  allowable  apogee  radius  be  larger  than 
the  Earth  radius.  Moreover,  the  second  coefficient  2(R  ■  u)  will  essentially  always  be  positive  because  observations 
have  to  be  taken  above  the  local  horizontal  plane  at  some  positive  local  elevation  angle. 

Also  as  before,  for  space-based  observing  stations,  it  is  possible  that  neither  of  these  circumstances  would  be  true: 
the  station’s  orbital  position  may  be  above  the  maximum  specified  apogee  radius,  or  observations  may  be  taken  at 
negative  local  elevation  angles,  or  both.  If  the  third  coefficient  in  (8)  is  positive,  that  is,  if  a^AX(l  ~l~  6max)2  <  R  ' 
R  ,  then  the  quadratic  will  have  either  no  positive  real  roots  or  two  positive  real  roots,  depending  on  the  sign  of  the 
second  coefficient.  This  is  the  possibility  just  mentioned  for  space-based  stations,  although  we  do  not  expect  this 
case  for  Earth-bound  stations.  If,  furthermore,  the  second  coefficient  in  (8)  is  positive,  that  is,  if  (R  ■  u)  >  0  ,  then 
we  have  no  positive  real  roots,  but  only  a  pair  of  negative  roots.  Because  the  quadratic  is  concave -up,  the  inequality 
(8)  is  satisfied  between  these  roots.  However,  since  we  require  a  priori  that  range  values  be  non-negative,  we  can 
discard  this  particular  observation  and  form  no  range  hypotheses  for  it. 


If  the  third  coefficient  in  (8)  is  positive,  but  the  second  coefficient  is  negative,  (R  ■  u)  <  0,  meaning  that  the 
observation  is  taken  at  negative  local  elevation  angle,  then  the  quadratic  will  have  two  positive  real  roots.  The 
quadratic  is  concave-up,  so  the  inequality  (8)  will  be  satisfied  between  these  two  roots.  In  this  case,  we  have  a  single 
finite  interval  of  range  over  which  range  hypotheses  will  satisfy  the  apogee  condition: 


p  >  -(R  ■  u) 


(R  ■  u)2  +  [a 


2 

MAX 


(1  +  eMAx)2  —  R  '  R] 


(11) 


p  <  -(R  ■  u)  + 


(R  ■  u)2  +  [a 


2 

MAX 


(1  +  £Max)2  —  R  '  R] 


(12) 


The  set  of  range  values  over  which  we  may  have  to  form  hypotheses  for  the  observation  in  question  is  given  by  the 
intersection  of  all  of  the  above  conditions,  both  perigee  conditions  and  apogee  conditions. 

The  above  conditions  are  bounds  on  the  possible  values  of  range,  which  can  be  computed  for  each  single 
observation.  The  fact  that  only  single  observations  are  involved  is  what  allows  us  to  find  explicit  bounds  for  each  of 
the  ranges  before  we  form  any  range  hypotheses.  However,  at  least  five  additional  restrictions  on  the  allowable 


values  of  range  can  be  deduced  from  relations  that  involve  both  of  the  ranges  presented  for  a  solution  to  Lambert’s 
problem.  Although  the  nonlinearities  in  these  relations  prevent  us  from  getting  explicit  inequalities  like  (5)  -  (7)  and 
(10)  -  (12),  nevertheless  we  can  formulate  additional  conditions  that  px  and  p2  must  satisfy.  Checking  these  extra 
conditions  for  each  range  pair  may  keep  us  from  having  to  produce  some  unnecessary  and  relatively  expensive 
Lambert  solutions. 

Using  the  vector  triangle  relation  for  each  of  the  two  lines  of  sight,  compute  the  unit  vector  n  normal  to  the 
candidate  orbital  plane: 

n  =  ±(rx  x  r2)/||r1  x  r2||  (13) 

Here  the  ambiguous  sign  is  resolved  as  “+”  for  “short- way”  trajectories  and  for  “long- way”  trajectories.  The 
inclination  is  given  unambiguously  by 

cos  /  =  n  ■  k  (14) 

The  inclination  of  the  candidate  orbit  lies  in  the  specified  interval  [/min  *  /max]  provided  that  cos  /  lies  in  the  interval 
[cos  /MAX  ,  cos  /MIN]  .  Hence  we  require  that 

cos /MAx  <  n  k  <  cos /Min  (15) 

In  a  similar  way,  we  use  the  unit  nodal  vector  to  obtain  conditions  that  the  range  pair  must  satisfy  if  the  candidate 
orbit  is  to  lie  within  a  specified  interval  of  right  ascension  of  the  ascending  node,  [HMin  >  ^max]  •  In  the  Earth- 
centered  inertial  frame,  we  have 

(k  x  n)/||k  x  n||  =  (cosfi,sinn,0)T  (16) 

so  that,  following  standard  logic  for  quadrant  resolution,  we  require 

HMin  ^  tan_1(sinn/cos  H)  <  HMax  (17) 

Of  course,  for  important  special  cases  like  near-GEO  orbits,  it  may  be  preferable  to  define  element-space  partitions 
in  terms  of  nonsingular  elements  such  as  p  =  sin(//2)cosH  and  q  =  sin(//2)sinH  .  No  special  difficulty 
attaches  to  working  in  terms  of  these  or  any  other  elements  related  to  the  orbit  plane.  If  any  range-pair  hypothesis 
(pi,p2)  does  not  satisfy  all  of  the  above  conditions,  then  that  pair  of  values  can  be  eliminated  from  further 
consideration  without  solving  Lambert’s  problem.  Note  that  it  is  the  pair  of  range  values  that  is  eliminated;  either 
range  value  by  itself  may  still  lead  to  an  acceptable  hypothesis  in  combination  with  some  other  range  value. 

Next,  we  can  use  three  special  solutions  of  Lambert’s  problem  to  restrict  the  ranges.  The  eccentricity  of  the  orbit  of 
least  possible  eccentricity  that  goes  through  a  given  pair  of  position  vectors  can  be  computed  solely  in  terms  of  those 
position  vectors.  Call  it  e0  .  Likewise,  the  semimajor  axis  of  the  orbit  of  least  possible  semimajor  axis  that  goes 
through  the  pair  of  positions  can  be  computed  solely  in  terms  of  the  position  vectors.  Call  it  a0  .  Hence,  for  each 
hypothesized  range  pair  (px  ,  p2)  ,  we  compute  the  corresponding  position  vectors  and  apply  the  following  logic: 

If  a0  >  aM ax  j  then  reject  the  hypothesis  pair  without  solving  Lambert’s  problem,  because  the 
geometry  is  guaranteed  to  produce  a  larger  semimajor  axis  than  we  have  specified. 

If  e0  >  £Max  ?  then  reject  the  hypothesis  pair  without  solving  Lambert’s  problem,  because  the 
geometry  is  guaranteed  to  produce  a  larger  eccentricity  than  we  have  specified. 

Of  course,  even  for  a  (px  ,p2)  hypothesis  that  passes  the  above  tests,  the  actual  solution  of  Lambert’s  problem  may 
still  turn  out  to  get  rejected  once  we  have  computed  the  elements  of  the  candidate  orbit.  The  reason  is  that  none  of 
the  conditions  on  range  derived  above  involves  the  minimum  allowable  eccentricity,  eMiN  •  This  fundamental  feature 
of  our  problem  raises  the  question  of  how  well  we  can  limit  the  generation  of  candidate  orbits  to  lie  within  the  given 
eccentricity  interval.  Let  us  assume  that  the  hypothetical  range  pair  is  not  rejected  by  the  above  criterion,  so  that 
e0  <  eMAX  •  Assume  also  that  all  of  the  range  bounds  and  other  conditions  that  depend  on  single  observations  have 
already  been  applied.  Then  we  know  that  the  Lambert  solution  for  a  pair  of  range  hypotheses  will  not  produce  an 
orbit  having  eccentricity  outside  the  interval  [e0  ,  eMAX]  .  If  eMiN  <  e0  ,  we  have  no  difficulty:  the  candidate  orbit 


will  have  an  eccentricity  within  the  given  interval  [eMIN  ,  eMAX]  •  However,  if  e0  <  eMIN  ,  then  the  eccentricity  of 
the  candidate  orbit  may  or  may  not  lie  within  the  specified  interval.  The  Lambert  solution  has  to  be  generated  and 
then  either  kept  if  the  eccentricity  is  at  least  as  large  as  eMIN  or  discarded  if  the  candidate  eccentricity  turns  out  to  be 
less  than  eMIN  .  This  represents  some  inefficiency  in  the  generation  of  candidate  orbits,  especially  if  those  same 
candidate  orbits  were  to  be  generated  in  the  processing  for  other  element-space  partitions.  The  extent  of  the  overall 
inefficiency  depends  on  the  dataset  and  the  actual  element-space  partitions  being  used,  so  we  cannot  draw  general 
conclusions.  It  would  be  helpful  at  this  point  to  have  reasonably  sharp  bounds  on  the  actual  eccentricity  in  the 
Lambert  problem  without  having  to  solve  the  whole  problem.  However,  lacking  that,  we  have  no  better  recourse 
than  to  generate  the  candidate  orbit.  Overall,  we  do  expect  to  be  able  to  reduce  the  number  of  Lambert  solutions  that 
have  to  be  generated,  compared  to  the  number  required  without  the  above  checks  involving  a0  and  e0. 


The  formulas  for  a0  and  e0  are  well  known: 


4a0  =  IM  +  llrJ  +  Hr.-rJ 


and 


= 


KIM-IM)! 
I|r2  -rj 


(18) 


Finally,  Euler’s  Theorem,  a  special  case  of  Lambert’s  Theorem,  expresses  the  time  of  flight  AtP  between  given 
position  vectors  on  a  parabolic  (zero-energy)  orbit: 


AtP 


4 

3 


(1  -  s  A3) 


(19) 


Here  the  quantity  5  is  a  signum  function:  s  =  +1  for  “short- way”  trajectories  and  s  =  —  1  for  “long-way” 
trajectories.  The  parameter  A  is  defined  in  terms  of  the  position  vectors: 


IM  +  IM-lte-rJ  (20) 

Hr,  ||  +  ||r2||  +  ||r2-r1||- 


Because,  for  given  position  vectors,  the  time  of  flight  in  Lambert’s  problem  is  a  monotonic  decreasing  function  of 
the  orbital  energy,  elliptic  (negative-energy)  orbits  will  always  have  a  time  of  flight  longer  than  the  parabolic  time, 
and  hyperbolic  (positive-energy)  orbits  will  always  have  a  time  of  flight  shorter  than  the  parabolic  time.  In  our  case, 
we  can  require  that  our  observation  pairs  and  range  hypotheses  always  produce  elliptic  orbits: 


^2  Tl  ^  Atp 


(21) 


Combinations  that  do  not  satisfy  this  condition  can  be  eliminated  without  generating  a  Lambert  solution.  Given  an 
observation  pair  ux  and  u2  ,  the  previous  formulas,  and  the  associated  logic,  can  be  used  to  decide  if  a  hypothetical 
pair  of  ranges  should  be  used  to  generate  a  Lambert  solution.  Of  course,  whatever  Lambert  solutions  are  generated 
should  be  verified  for  compliance  with  the  specified  interval  of  eccentricity,  because  none  of  the  conditions  on  range 
derived  so  far  depends  on  the  value  of  the  minimum  allowable  eccentricity  eMIN  .  By  construction,  we  have 
guaranteed  compliance  with  the  intervals  of  semimajor  axis,  inclination  and  right  ascension  of  the  ascending  node. 


Figure  2:  (a)  Raw  unidentified  observations  (b)  Observations  that  have  been  randomly  connected 
together 


3.  APPLYING  RANGE  BOUNDARIES  TO 
RANGE  ESTIMATION 

Now  that  we  have  defined  a  set  of  rigorous  boundaries  for 
candidate  range  estimates,  we  can  combine  these  boundaries 
with  traditional  range  estimation  procedures  for  angles-only 
IOD.  First,  blindly  generate  all  possible  combinations  of 
observations  (see  Figure  2).  Next,  using  the  method  of  Gauss, 
estimate  ranges  for  each  of  the  observations  within  a  given 
combination  using  the  techniques  outlined  in  [5]  and  [9].  It  is 
not  necessary  to  finish  the  IOD  process.  At  this  stage  we  only 
calculate  the  range  estimates  for  each  observation  in  a  draw. 
Filter  these  range  estimates  for  obviously  invalid  solutions 
before  applying  the  orbit  element  boundaries  to  further 
eliminate  infeasible  solutions.  The  last  step  is  to  check  range 
estimates  for  consistency  and  persistence  across  combinations 
returning  potentially  valid  solutions. 


■ 

Obs. 

Indices 

Pi 

P2 

P3 

1 

(1,10,19) 

412 

457 

435 

PASS 

2 

(1,10,28) 

420 

475 

490 

PASS 

3 

(1,17,28) 

364 

621 

400 

FAIL 

m 

(2,11,29) 

537 

545 

523 

PASS 

m+1 

(2,21,35) 

678 

323 

721 

FAIL 

Table  2:  Estimates  of  range  for  various  random  draws  of 
observations. 

Here  we  introduce  the  concept  of  consistency  and  persistence. 
These  concepts  hinge  upon  the  notion  that  all  IOD  range 
estimates  will  be  inaccurate  and  unrefined  but  generated  using 
a  common  method.  Using  consistency  and  persistence,  we 
will  assign  the  random  draws  of  observations  into  one  of  three 
categories:  Not  Likely ,  Likely ,  and  More  Likely.  Obviously, 
failed  cases  automatically  go  into  the  Not  Likely  category. 


Observation  ID  Numbers 

1 

10 

19 

28 

2 

11 

18 

20 

29 

3 

12 

19 

21 

30 

32 

4 

13 

21 

22 

31 

33 

5 

14 

23 

32 

6 

15 

24 

33 

7 

16 

25 

34 

8 

17 

26 

35 

9 

18 

27 

36 

10 

1 

19 

28 

11 

2 

20 

27 

29 

32 

12 

3 

21 

30 

13 

4 

22 

31 

14 

5 

23 

32 

15 

6 

24 

33 

16 

7 

25 

34 

17 

8 

19 

26 

32 

35 

18 

2 

9 

22 

27 

36 

19 

1 

3 

10 

17 

28  32 

20 

2 

11 

29 

21 

3 

4 

12 

30 

33 

22 

2 

4 

13 

18 

31 

23 

5 

14 

32 

24 

6 

15 

33 

25 

7 

16 

34 

26 

8 

17 

35 

27 

9 

11 

18 

32 

36 

28 

1 

10 

19 

29 

2 

11 

20 

30 

3 

12 

21 

31 

4 

13 

22 

32 

3 

5 

11 

14 

17  19  23  27 

33 

4 

6 

15 

21 

24 

34 

7 

16 

25 

35 

8 

17 

26 

36 

9 

18 

27 

Table  1:  Persistence  table  where  all  valid  associations 
between  individual  observations  can  be  seen. 


There  are  two  types  of  consistency.  The  first  is  consistency  of  range  estimates  within  a  random  draw.  For  cases 
where  the  normalized  time  between  observations  is  small,  then  physics  dictates  that  the  range  estimates  across  the 
three  observations  would  be  roughly  similar.  The  second  is  consistency  of  range  estimates  for  a  specific  observation 
across  all  random  draws.  Many  bad  cases  will  be  caught  by  the  orbit  element  boundaries  but  one  can  also  establish  a 
mean  and  variance  for  the  estimate  of  range  associated  with  each  observation.  These  two  types  of  consistency  can  be 
used  to  distinguish  outlier  range  estimates  if  there  is  trouble  discerning  which  of  several  valid  random  draws  should 
be  flagged  as  more  or  less  likely. 


Consistent  solutions  can  coalesce  into  persistent  solutions.  A  persistent  solution  occurs  when  a  particular 
observation  is  repeatedly  associated  with  another  observation  in  the  list  of  valid  random  draws.  Table  1  shows  a  list 
of  hypothetical  range  estimates  for  different  sets  of  observation  combinations.  Individual  cases  are  first  screened  for 
validity  based  upon  the  orbit  element  boundaries  described  previously  and  are  marked  pass  or  fail.  Consider  the 


Ranking 


observation  triplets  (1,10,19)  and  (10,19,28)  that  have  both  passed  the  orbit 
element  boundary  tests  meaning  that  these  are  Likely  observation  pairings.  We  can 
check  the  range  estimates  for  observation  10  and  19  for  consistency  if  desired  to 
further  distinguish  observation  groupings.  However,  persistence  says  that  if 
(1,10,19)  as  well  as  (10,19,28)  are  likely  valid,  then  (1,10,19,28)  must  also  be  a 
Likely  grouping  of  observations. 

Following  this  process  for  all  valid  observation  draws,  we  can  generate  Table  2 
where  our  initial  valid  groupings  of  three  observations  have  been  strung  together 
into  groups  of  four  or  more  observations.  Looking  at  row  one,  we  see  that 
observation  1  associated  only  with  10,  19,  and  28.  Furthermore,  row  10  indicates 
that  observation  10  associated  only  with  1,19,  28.  Likewise,  row  28  indicates  that 
observation  28  associated  only  with  1,  10,  and  19.  However,  observation  19  on 
row  19  associated  with  1,3,  10,  17,  28,  32.  Using  persistence,  we  can  say  that  the 
association  between  1,  10,  19,  and  28  is  More  Likely  because  3  out  of  the  4 
observations  associated  only  with  each  other  and  the  fourth  observation  associated 
with  a  few  additional  observations.  The  question  then  is  should  we  toss  out  the 
candidate  grouping  of  (19,  1,3,  10,  17,  28,  32)?  To  be  conservative  one  can  flag 
this  case  as  Likely  and  maintain  it  in  a  multiple  hypothesis  setting.  One  could  also 
look  at  the  consistency  checks  for  each  range  estimate  to  make  a  further 
discrimination  of  likelihood.  However,  further  inspection  of  Table  2  reveals  that 
observations  3,  17,  and  32  are  more  strongly  associated  with  other  observations. 
Therefore,  one  could  safely  demote  this  candidate  grouping  in  the  rankings. 

Proceeding  through  each  row  of  Table  2,  we  generate  a  ranking  of  hypothetical  observation  groupings  in  Table  3. 
Cases  where  certain  groupings  of  observations  associated  only  with  themselves  are  flagged  as  More  Likely.  The  one- 
off  cases  of  extended  groupings  can  be  maintained  as  Likely  or  demoted  in  the  rankings.  In  any  case,  one  should 
always  maintain  the  original  set  of  valid  three-observation  groupings  as  Likely  in  addition  to  the  built-up  groupings. 
All  of  the  final  ranked  observation  groupings  can  now  be  used  to  generate  candidate  orbits  in  a  multiple  hypothesis 
sense  to  examine  refined  orbit  solutions.  Essentially  this  use  of  persistence  and  consistency  allows  us  to  perform  a 
“poor  man’s  version”  of  the  classical  assignment  problem.  Future  work  will  attempt  to  employ  the  formal 
methodologies  of  the  classical  assignment  problem  with  the  range  estimate  consistency  checks  as  a  cost  function. 
Our  hope  would  be  to  determine  a  more  nearly  optimal  set  of  observation  assignments. 

4.  A  MONTE  CARLO  SIMULATION 

Using  a  simple  two-body  orbit  propagator,  angles-only  observations  of  Intelsat-8  were  simulated  as  viewed  from 
PanSTARRS  on  top  of  Mt.  Haleakala  (203.74°  longitude,  20.71°  latitude,  3.07  km  elevation  about  the  WGS-95 
geoid).  Observations  were  obtained  two  hours  apart  over  a  24  hour  period  for  a  total  of  13  observations.  One 
observation  was  selectively  “bumped”  with  a  random  Gaussian  N(  120, 120)  arc-second  error.  Using  the  above 
outlined  approached,  each  of  the  13  observations  was  thrown  into  a  “bucket”  and  drawn  blindly  in  groupings  of  4 
observations  at  a  time  leading  to  a  possible  715  combinations.  Of  these  715  combinations,  220  combinations 
contained  the  “bumped”  or  “bad”  observation.  The  question  to  be  answered  is  “Can  we  prune  out  the  hypotheses 
with  the  bad  observation  a  prioriT ’  Using  the  range  estimation  procedure  detailed  above,  the  range  estimates  were 
computed  for  each  of  the  715  combinations  of  4  observations.  The  physics-based  boundaries  were  tailored  for  a 
GEO  object  and  set  at  semi-major  axis  in  the  interval  [37500,  45000]  km,  eccentricity  in  [0,  0.075],  and  inclination 
in  [-12°,  12°].  The  range  estimates  could  fail  for  a  variety  of  reasons  including  poor  geometry  and  negative  range 
estimates.  Only  hypotheses  that  returned  a  valid  solution  were  subjected  to  the  orbit  element  boundary  filters.  This 
process  was  repeated  in  a  Monte  Carlo  fashion  for  1000  trials  with  a  specified  observation  having  random  error 
applied  to  it.  The  results  of  this  simulation  were  that  out  of  the  715,000  total  hypotheses  over  1000  trials,  only  2968 
cases  passed  through  the  various  validation  procedures  and  resulted  in  a  false  association.  This  means  that  we  had  a 
0.42%  acceptance  rate  for  the  bad  observation. 


Obs 


More  Likely  (1,10,19,28) 

More  Likely  (2,11,20,29) 

More  Likely  (3,12,21,30) 

Likely  (1,10,19) 

Likely  (2,11,20) 

Not  Likely  (1,10,33) 

Not  Likely  (2,21,35) 

Not  Likely  (5,17,24) 

Table  3:  Ranking  of  the 
persistent  and  consistent 
solutions. 


5.  A  GEO  CLUSTER  SIMULATION 


Now  we  consider  the  case  of  a  GEO  cluster  of  9  satellites  separated  by  0.25°  in  right  ascension  of  the  ascending 
node.  For  each  of  the  nine  objects,  four  simulated  observations  2  hours  apart  were  generated  with  2  arc  seconds  of 
error  applied.  In  total  there  were  36  observations  that  were  drawn  blindly  3  at  a  time  for  a  total  of  7140  possible 
combinations.  The  question  now  is  “Can  we  rank  hypotheses  in  order  of  most  likely  to  least  likely  a  prioriT ’  The 
result  of  this  simulation  is  that,  after  computing  the  range  estimates  and  applying  the  physics  based  boundaries  on 
orbit  elements,  the  algorithm  returned  valid  range  estimates  for  all  of  the  cases  that  we  knew  to  be  correct  plus  a 
handful  of  cases  that  were  incorrect.  Using  consistency  and  persistence  to  string  together  observations  into  longer 
arcs  of  data,  the  algorithm  correctly  built  up  the  3 -observation  groupings  into  the  correct  4-observation  groupings 
and  only  these  cases  were  marked  as  highly  likely.  None  of  the  incorrect  3 -observation  pairings  from  the  previous 
step  were  included  in  the  longer  arcs  of  data  and  these  were  all  flagged  as  less  likely  due  to  their  unique  occurrences. 

6.  CONCLUSION 

Our  results  show  that  the  possible  values  of  both  range  and  range  rate  can  be  limited  a  priori  for  each  line-of-sight 
observation  to  finite  intervals  corresponding  to  a  specified  partition  of  the  element  space.  The  endpoints  of  the 
intervals  are  given  explicitly  in  terms  of  the  angle-based  observations,  station  position  and  station  velocity,  and  can 
be  computed  independently  for  each  observation.  Additional  conditions  based  on  the  orientation  of  the  orbital  plane 
and  special  solutions  of  Lambert’s  problem,  which  must  be  satisfied  by  range  values  for  pairs  of  observations,  can 
be  used  to  further  reduce  the  number  of  Lambert  solutions  needed  for  the  initial  orbit  determinations.  All  the 
formulas  derived  here  apply  uniformly  to  Earth-bound  and  space-based  observing  stations.  We  also  present  explicit 
conditions  identifying  when  a  given  observation  does  not  correspond  to  any  possible  orbit  within  the  specified 
element-space  partition.  Such  observations  can  be  discarded  before  any  data  association  hypotheses  or  Lambert 
solutions  are  produced.  We  have  demonstrated  the  effectiveness  of  these  techniques  through  a  single  satellite  Monte 
Carlo  simulation  and  a  GEO  cluster  simulation.  We  have  found  that  it  is  possible  to  rank  observation  hypotheses  in 
rough  order  of  likelihood  a  priori.  We  have  found  this  method  is  also  effective  for  culling  out  bad  observations. 
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